#opening csv
soils <- read_csv(here("data", "shift_soil_data.csv")) %>% # read in soil csv
mutate(date = lubridate::mdy(date), # change to date class
week = lubridate::week(date)) %>% # create a week column
select(soil_id:transect, meter_location, electro_cond_mS_per_cm, season, latitude, longitude, landcover)%>% #select certain column to work with
mutate(paired_flight = lubridate::ymd(paired_flight)) %>%
drop_na() %>% #drop any NA rows
st_as_sf(coords = c("longitude", "latitude")) %>%
st_set_crs(value = 4326) %>%
st_transform(crs = 32611)
# Open spatial Data
soils_sf <- read_sf(here("data", "soil_w_elevation.shp")) %>% # read in sf file
mutate(date = lubridate::ymd(date), #change date to date class
week = lubridate::week(date)) %>% #create week
select(electro_co, date, transect, meter_loca) %>% # select certain columns
filter(date != "2022-09-12") %>%
mutate(electro_co = as.numeric(electro_co)) %>% #change to numeric class
drop_na() #drop NAs
bbox <- read_sf(here("data", "bbox.shp")) #read in bounding box file
boundary <- read_sf(here("data", "marsh_boundary.shp")) %>%
st_transform(crs = 32611)
dem <- read_stars(here("data", "dem.tif")) %>% # read in dem
st_crop(y = st_bbox(bbox)) %>% #crop to bounding box
st_warp(cellsize = 4.8, crs = 32611)
mllw <- dem + 0.042
mllw_all <- c(mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
mllw_extract <- mllw %>% # call dem
st_extract(soils$geometry) # extract elevation to point layer
soils <- cbind(soils, mllw_extract$dem.tif) %>% #bind extracted values to the soil data
rename("elevation" = "mllw_extract.dem.tif") %>%
st_transform(crs = 32611)
savi_1 <- read_stars(here("imagery", "subset_images", "rfl_2022_02_24.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_2 <- read_stars(here("imagery", "subset_images", "rfl_2022_02_28.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_3 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_08.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_4 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_16.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_5 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_22.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_6 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_05.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_7 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_12.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_8 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_20.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_9 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_29.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_10 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_03.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_11 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_11.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_13 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_17.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_14 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_29.tif")) %>%
st_set_dimensions(3, values = seq(1:425)) %>%
split(3) %>%
mutate(savi = ((((X85+X86)/2)-((X59+X60)/2))/(((X85+X86)/2)+((X59+X60)/2)+0.5)) * 1.5) %>%
select(savi)
savi_all <- c(savi_1, savi_2, savi_3, savi_4, savi_5, savi_6, savi_7, savi_8, savi_9, savi_10, savi_11, savi_13, savi_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
ndvi_1 <- read_stars(here("imagery", "NDVI", "2022_02_24_ndvi.tif")) %>%
setNames("ndvi")
ndvi_2 <- read_stars(here("imagery", "NDVI", "2022_02_28_ndvi.tif")) %>%
setNames("ndvi")
ndvi_3 <- read_stars(here("imagery", "NDVI", "2022_03_08_ndvi.tif")) %>%
setNames("ndvi")
ndvi_4 <- read_stars(here("imagery", "NDVI", "2022_03_16_ndvi.tif")) %>%
setNames("ndvi")
ndvi_5 <- read_stars(here("imagery", "NDVI", "2022_03_22_ndvi.tif")) %>%
setNames("ndvi")
ndvi_6 <- read_stars(here("imagery", "NDVI", "2022_04_05_ndvi.tif")) %>%
setNames("ndvi")
ndvi_7 <- read_stars(here("imagery", "NDVI", "2022_04_12_ndvi.tif")) %>%
setNames("ndvi")
ndvi_8 <- read_stars(here("imagery", "NDVI", "2022_04_20_ndvi.tif")) %>%
setNames("ndvi")
ndvi_9 <- read_stars(here("imagery", "NDVI", "2022_04_29_ndvi.tif")) %>%
setNames("ndvi")
ndvi_10 <- read_stars(here("imagery", "NDVI", "2022_05_03_ndvi.tif")) %>%
setNames("ndvi")
ndvi_11 <- read_stars(here("imagery", "NDVI", "2022_05_11_ndvi.tif")) %>%
setNames("ndvi")
#ndvi_12 <- read_stars(here("imagery", "NDVI", "2022_05_12_ndvi.tif")) %>%
# setNames("ndvi")
ndvi_13 <- read_stars(here("imagery", "NDVI", "2022_05_17_ndvi.tif")) %>%
setNames("ndvi")
ndvi_14 <- read_stars(here("imagery", "NDVI", "2022_05_29_ndvi.tif")) %>%
setNames("ndvi")
ndvi_all <- c(ndvi_1, ndvi_2, ndvi_3, ndvi_4, ndvi_5, ndvi_6, ndvi_7, ndvi_8, ndvi_9, ndvi_10, ndvi_11, ndvi_13, ndvi_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
vssi_1 <- read_stars(here("imagery", "VSSI", "02_24_vssi.tif")) %>%
setNames("vssi")
vssi_2 <- read_stars(here("imagery", "VSSI", "02_28_vssi.tif")) %>%
setNames("vssi")
vssi_3 <- read_stars(here("imagery", "VSSI", "03_08_vssi.tif")) %>%
setNames("vssi")
vssi_4 <- read_stars(here("imagery", "VSSI", "03_16_vssi.tif")) %>%
setNames("vssi")
vssi_5 <- read_stars(here("imagery", "VSSI", "03_22_vssi.tif")) %>%
setNames("vssi")
vssi_6 <- read_stars(here("imagery", "VSSI", "04_05_vssi.tif")) %>%
setNames("vssi")
vssi_7 <- read_stars(here("imagery", "VSSI", "04_12_vssi.tif")) %>%
setNames("vssi")
vssi_8 <- read_stars(here("imagery", "VSSI", "04_20_vssi.tif")) %>%
setNames("vssi")
vssi_9 <- read_stars(here("imagery", "VSSI", "04_29_vssi.tif")) %>%
setNames("vssi")
vssi_10 <- read_stars(here("imagery", "VSSI", "05_03_vssi.tif")) %>%
setNames("vssi")
vssi_11 <- read_stars(here("imagery", "VSSI", "05_11_vssi.tif")) %>%
setNames("vssi")
#vssi_12 <- read_stars(here("imagery", "VSSI", "05_12_vssi.tif")) %>%
# setNames("vssi")
vssi_13 <- read_stars(here("imagery", "VSSI", "05_17_vssi.tif")) %>%
setNames("vssi")
vssi_14 <- read_stars(here("imagery", "VSSI", "05_29_vssi")) %>%
setNames("vssi")
vssi_all <- c(vssi_1, vssi_2, vssi_3, vssi_4, vssi_5, vssi_6, vssi_7, vssi_8, vssi_9, vssi_10, vssi_11, vssi_13, vssi_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
mari_1 <- read_stars(here("imagery", "mARI", "02_24_mari_avg.tif")) %>%
setNames("mari")
mari_2 <- read_stars(here("imagery", "mARI", "02_28_mari_avg.tif")) %>%
setNames("mari")
mari_3 <- read_stars(here("imagery", "mARI", "03_08_mari_avg.tif")) %>%
setNames("mari")
mari_4 <- read_stars(here("imagery", "mARI", "03_16_mari_avg.tif")) %>%
setNames("mari")
mari_5 <- read_stars(here("imagery", "mARI", "03_22_mari_avg.tif")) %>%
setNames("mari")
mari_6 <- read_stars(here("imagery", "mARI", "04_05_mari_avg.tif")) %>%
setNames("mari")
mari_7 <- read_stars(here("imagery", "mARI", "04_12_mari_avg.tif")) %>%
setNames("mari")
mari_8 <- read_stars(here("imagery", "mARI", "04_20_mari_avg.tif")) %>%
setNames("mari")
mari_9 <- read_stars(here("imagery", "mARI", "04_29_mari_avg.tif")) %>%
setNames("mari")
mari_10 <- read_stars(here("imagery", "mARI", "05_03_mari_avg.tif")) %>%
setNames("mari")
mari_11 <- read_stars(here("imagery", "mARI", "05_11_mari_avg.tif")) %>%
setNames("mari")
#mari_12 <- read_stars(here("imagery", "mARI", "05_12_mari_avg.tif")) %>%
# setNames("mari")
mari_13 <- read_stars(here("imagery", "mARI", "05_17_mari_avg.tif")) %>%
setNames("mari")
mari_14 <- read_stars(here("imagery", "mARI", "05_29_mari_avg.tif")) %>%
setNames("mari")
mari_all <- c(mari_1, mari_2, mari_3, mari_4, mari_5, mari_6, mari_7, mari_8, mari_9, mari_10, mari_11, mari_13, mari_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
crsi_1 <- read_stars(here("imagery", "CRSI", "02_24_crsi.tif")) %>%
setNames("crsi")
crsi_2 <- read_stars(here("imagery", "CRSI", "02_28_crsi.tif")) %>%
setNames("crsi")
crsi_3 <- read_stars(here("imagery", "CRSI", "03_08_crsi.tif")) %>%
setNames("crsi")
crsi_4 <- read_stars(here("imagery", "CRSI", "03_16_crsi.tif")) %>%
setNames("crsi")
crsi_5 <- read_stars(here("imagery", "CRSI", "03_22_crsi.tif")) %>%
setNames("crsi")
crsi_6 <- read_stars(here("imagery", "CRSI", "04_05_crsi.tif")) %>%
setNames("crsi")
crsi_7 <- read_stars(here("imagery", "CRSI", "04_12_crsi.tif")) %>%
setNames("crsi")
crsi_8 <- read_stars(here("imagery", "CRSI", "04_20_crsi.tif")) %>%
setNames("crsi")
crsi_9 <- read_stars(here("imagery", "CRSI", "04_29_crsi.tif")) %>%
setNames("crsi")
crsi_10 <- read_stars(here("imagery", "CRSI", "05_03_crsi.tif")) %>%
setNames("crsi")
crsi_11 <- read_stars(here("imagery", "CRSI", "05_11_crsi.tif")) %>%
setNames("crsi")
#crsi_12 <- read_stars(here("imagery", "CRSI", "05_12_crsi.tif")) %>%
# setNames("crsi")
crsi_13 <- read_stars(here("imagery", "CRSI", "05_17_crsi.tif")) %>%
setNames("crsi")
crsi_14 <- read_stars(here("imagery", "CRSI", "05_29_crsi.tif")) %>%
setNames("crsi")
crsi_all <- c(crsi_1, crsi_2, crsi_3, crsi_4, crsi_5, crsi_6, crsi_7, crsi_8, crsi_9, crsi_10, crsi_11, crsi_13, crsi_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
frac_1 <- read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_02_24_shadeNorm")) %>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_2 <-read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_02_28_shadeNorm")) %>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_3 <-read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_03_08_shadeNorm"))%>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_4 <-read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_03_16_shadeNorm"))%>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_5 <-read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_03_22_shadeNorm")) %>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_6 <- read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_04_05_shadeNorm"))%>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_7 <- read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_04_12_shadeNorm"))%>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_8 <- read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_04_20_shadeNorm"))%>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_9 <- read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_04_29_shadeNorm"))%>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_10 <- read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_05_03_shadeNorm"))%>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_11 <- read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_05_11_shadeNorm"))%>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_13 <- read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_05_17_shadeNorm"))%>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_14 <- read_stars(here("imagery", "fractional_cover", "shade_norm", "rfl_2022_05_29_shadeNorm"))%>%
filter(band %in% c("ShadeNorm_ NPV_fraction", "ShadeNorm_ GREEN VEG_fraction", "ShadeNorm_ BARE SOIL_fraction")) %>%
split(3) %>%
setNames(c("Soil_fraction", "GV_fraction", "NPV_fraction" ))
frac_all <- c(frac_1, frac_2, frac_3, frac_4, frac_5, frac_6, frac_7, frac_8, frac_9, frac_10, frac_11, frac_13, frac_14, along = 3) %>%
st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")
npv_frac_all <- frac_all %>%
select("NPV_fraction")
gv_frac_all <- frac_all %>%
select("GV_fraction")
soil_frac_all <- frac_all %>%
select("Soil_fraction")
mllw_all <- mllw_all %>%
st_warp(dest = crsi_all) %>%
setNames("elevation")
savi_all <- savi_all %>%
st_warp(dest = crsi_all) %>%
st_as_stars()
soils_subset <- soils %>%
filter(paired_flight %in% lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")))
soils_extract_2 <- data.frame()
for (i in unique(lubridate::week(soils_subset$paired_flight))){
ndvi_week <- ndvi_all %>%
filter(lubridate::week(date) == i)
mari_week <- mari_all %>%
filter(lubridate::week(date) == i)
vssi_week <- vssi_all %>%
filter(lubridate::week(date) == i)
crsi_week <- crsi_all %>%
filter(lubridate::week(date) == i)
frac_week <- frac_all %>%
filter(lubridate::week(date) == i) %>%
st_as_sf()
savi_week <- savi_all %>%
filter(lubridate::week(date) == i)
soil_week <- soils %>%
mutate(week = lubridate::week(paired_flight)) %>%
filter(week == i)
ndvi_extract <- ndvi_week %>%
st_extract(soil_week)
mari_extract <- mari_week %>%
st_extract(soil_week)
vssi_extract <- vssi_week %>%
st_extract(soil_week)
crsi_extract <- crsi_week %>%
st_extract(soil_week)
frac_extract <- soil_week %>%
st_join(frac_week) %>%
select("Soil_fraction":"NPV_fraction") %>%
st_drop_geometry()
savi_extract <- savi_week %>%
st_extract(soil_week)
soils_week <- bind_cols(soil_week, ndvi_extract$ndvi, mari_extract$mari, vssi_extract$vssi, crsi_extract$crsi, savi_extract$savi, frac_extract) %>%
rename("ndvi" = "...12",
"mari" = "...13",
"vssi" = "...14",
"crsi" = "...15",
"savi" = "...16")
soils_extract_2 <- rbind(soils_extract_2, soils_week)
}
#write.csv(soils_extract, file = here("data", "extracted_soils.csv"))
set.seed(1234)
fit_all_savi <- train(electro_cond_mS_per_cm ~ elevation + mari + vssi + savi,
data = soils_extract_2,
method = "rf",
trControl = trainControl(method = "cv", number = 10),
importance = TRUE)
# using caret
fit_all_savi_imp <- as.data.frame(fit_all_savi$finalModel$importance)
fit_all_savi_imp_scaled <- predict(preProcess(fit_all_savi_imp, method = c("range")), fit_all_savi_imp) %>%
rownames_to_column(var = "variable")
ggplot(fit_all_savi_imp_scaled) +
geom_point(aes(x = IncNodePurity, y = variable), color = "#E69512")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = IncNodePurity), color = "#E69512", size = 1.2)+
labs(y = "Variables",#labels
x = "Importance (Scaled)")+
ggtitle("Increase in Node Purity")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 24, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 16),
axis.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
ggplot(fit_all_savi_imp) +
geom_point(aes(x = fit_all_savi_imp$`%IncMSE`, y = row.names(fit_all_savi_imp)), color ="#3A5D3D")+
geom_segment(aes(y = row.names(fit_all_savi_imp), x = 0, yend = row.names(fit_all_savi_imp), xend = fit_all_savi_imp$`%IncMSE`), color = "#3A5D3D", linetype = "dashed")+
labs(y = "Variables",#labels
x = "Importance (%)")+
ggtitle("% Increase in MSE")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 24, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 16),
axis.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "all_imp.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
ggplot() +#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted") +
geom_point(aes(y = fit_all_savi$finalModel$y, x = fit_all_savi$finalModel$predicted), color = "#E69512") +#data to make points from
geom_smooth(method = "lm", aes(y = fit_all_savi$finalModel$y, x = fit_all_savi$finalModel$predicted), color = "#3A5D3D") +#smooth line data
scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 14)) +
labs(y = "Actual Electrical Conductivity (mS/cm)", #labels
x = "Predicted Electrical Conductivity (mS/cm)") +
ggtitle("Elevation and Indices") +#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 24, hjust = 0.5), #plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 16),
axis.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "all_plot_caret.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
all_corr <- cor(y= fit_all_savi$finalModel$y, x = fit_all_savi$finalModel$predicted)
all_corr
## [1] 0.8502606
lm(fit_all_savi$finalModel$y ~ fit_all_savi$finalModel$predicted)
##
## Call:
## lm(formula = fit_all_savi$finalModel$y ~ fit_all_savi$finalModel$predicted)
##
## Coefficients:
## (Intercept) fit_all_savi$finalModel$predicted
## 0.04469 1.00473
max(fit_all_savi$finalModel$predicted)
## [1] 10.71859
min(fit_all_savi$finalModel$predicted)
## [1] 0.5540848
raster_all <- c(crsi_all, ndvi_all, vssi_all, mari_all, mllw_all, savi_all)
all_pred_savi <- predict(raster_all, fit_all_savi, drop_dimensions = F, na.action = na.omit)
#write_stars(all_pred_savi, dsn = here("data", "prediction_savi.tif"))
all_pred_crop <- all_pred_savi %>%
st_crop(y = boundary, crop = TRUE, as_points = FALSE) %>%
filter(date == "2022-02-24")
ggplot()+
geom_stars(data = all_pred_crop, aes(x = x, y = y, fill = prediction), na.action = na.omit) +
scale_fill_gradientn(colors = rev(cal_palette(name = "kelp1", n = 4, type = "continuous")))+
coord_sf(crs = 32116)+
labs(y = "Latitude",#labels
x = "Longitude",
fill = "Electrical
Conductivity
(mS/cm)")+
ggtitle("February 24, 2022")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"))
set.seed(1234)
fit_all_nc <- train(electro_cond_mS_per_cm ~ elevation + mari + vssi + ndvi,
data = soils_extract_2,
method = "rf",
trControl = trainControl(method = "cv", number = 10),
importance = TRUE)
# using caret
fit_all_nc_imp <- as.data.frame(fit_all_nc$finalModel$importance)
fit_all_nc_imp_scaled <- predict(preProcess(fit_all_nc_imp, method = c("range")), fit_all_nc_imp) %>%
rownames_to_column(var = "variable")
ggplot(fit_all_nc_imp_scaled) +
geom_point(aes(x = IncNodePurity, y = variable), color = "#E69512")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = IncNodePurity), color = "#E69512", size = 1.2)+
geom_point(aes(x = fit_all_nc_imp_scaled$`%IncMSE`, y = variable), color ="#3A5D3D")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = fit_all_nc_imp_scaled$`%IncMSE`), color = "#3A5D3D", linetype = "dashed")+
labs(y = "Variables",#labels
x = "Importance (Scaled)")+
ggtitle("Random Forest Variable Importance")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 24, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 16),
axis.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "all_imp.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(aes(y = fit_all_nc$finalModel$y, x = fit_all_nc$finalModel$predicted), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", aes(y= fit_all_nc$finalModel$y, x = fit_all_nc$finalModel$predicted), color = "#3A5D3D")+#smooth line data
scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Elevation and Spectral Indices")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 24, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 16),
axis.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
ggsave(filename = "all_plot_caret.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
all_corr <- cor(y= fit_all_nc$finalModel$y, x = fit_all_nc$finalModel$predicted)
all_corr
## [1] 0.8431511
lm(fit_all_nc$finalModel$y ~ fit_all_nc$finalModel$predicted)
##
## Call:
## lm(formula = fit_all_nc$finalModel$y ~ fit_all_nc$finalModel$predicted)
##
## Coefficients:
## (Intercept) fit_all_nc$finalModel$predicted
## 0.01662 1.00271
max(fit_all_nc$finalModel$predicted)
## [1] 10.64588
min(fit_all_nc$finalModel$predicted)
## [1] 0.5456866
all_pred_nc <- predict(raster_all, fit_all_nc, drop_dimensions = F, na.action = na.omit)
plot(all_pred_nc)
all_pred_crop <- all_pred_nc %>%
st_crop(y = boundary, crop = TRUE, as_points = FALSE) %>%
filter(date == "2022-02-24")
#write_stars(all_pred_crop, dsn = here("data", "prediction_ndvi.tif"))
ggplot()+
geom_stars(data = all_pred_crop, aes(x = x, y = y, fill = prediction), na.action = na.omit) +
scale_fill_gradientn(colors = rev(cal_palette(name = "kelp1", n = 4, type = "continuous")))+
coord_sf(crs = 32116)+
labs(y = "Latitude",#labels
x = "Longitude",
fill = "Electrical
Conductivity
(mS/cm)")+
ggtitle("February 24, 2022")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"))
t.test(x = fit_all_nc$finalModel$predicted, y = fit_all_savi$finalModel$predicted)
##
## Welch Two Sample t-test
##
## data: fit_all_nc$finalModel$predicted and fit_all_savi$finalModel$predicted
## t = 0.1192, df = 293.99, p-value = 0.9052
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.5595717 0.6317279
## sample estimates:
## mean of x mean of y
## 4.045458 4.009380
effsize::cohen.d(fit_all_nc$finalModel$predicted, fit_all_savi$finalModel$predicted)
##
## Cohen's d
##
## d estimate: 0.01385722 (negligible)
## 95 percent confidence interval:
## lower upper
## -0.2149285 0.2426430
model <- lm(electro_cond_mS_per_cm ~ elevation + mari + vssi + savi, data = soils_extract_2)
car::vif(model)
## elevation mari vssi savi
## 1.227227 2.417652 1.335792 2.207990
model <- lm(electro_cond_mS_per_cm ~ elevation + mari + vssi + savi + GV_fraction, data = soils_extract_2)
car::vif(model)
## elevation mari vssi savi GV_fraction
## 1.291167 2.489293 1.672837 2.479424 2.390923
No real change in GV fraction, Soil is still high, perhaps it was soil fractions in actuality
model <- lm(electro_cond_mS_per_cm ~ elevation + mari + vssi + ndvi + NPV_fraction + GV_fraction, data = soils_extract_2)
car::vif(model)
## elevation mari vssi ndvi NPV_fraction GV_fraction
## 1.332441 3.325679 1.505803 3.932088 1.209710 2.684775
GV dropped but MARI and NDVI increased, lets remove GV and/or ndvi
# remove ndvi
model <- lm(electro_cond_mS_per_cm ~ elevation + mari + vssi + GV_fraction, data = soils_extract_2)
car::vif(model)
## elevation mari vssi GV_fraction
## 1.261273 1.617618 1.467672 2.129178
# Add Soil remove GV
model <- lm(electro_cond_mS_per_cm ~ elevation + mari + vssi + Soil_fraction, data = soils_extract_2)
car::vif(model)
## elevation mari vssi Soil_fraction
## 1.231796 1.611226 1.757564 2.515505
removing ndvi and leaving GV fractions leaves acceptable collinearity
set.seed(1234)
fit_all_frac <- train(electro_cond_mS_per_cm ~ elevation + mari + vssi + Soil_fraction,
data = soils_extract_2,
method = "rf",
trControl = trainControl(method = "cv", number = 10),
importance = TRUE)
# using caret
fit_all_frac_imp <- as.data.frame(fit_all_frac$finalModel$importance)
fit_all_frac_imp_scaled <- predict(preProcess(fit_all_frac_imp, method = c("range")), fit_all_frac_imp) %>%
rownames_to_column(var = "variable")
ggplot(fit_all_frac_imp_scaled) +
geom_point(aes(x = IncNodePurity, y = variable), color = "#E69512")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = IncNodePurity), color = "#E69512", size = 1.2)+
geom_point(aes(x = fit_all_frac_imp_scaled$`%IncMSE`, y = variable), color ="#3A5D3D")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = fit_all_frac_imp_scaled$`%IncMSE`), color = "#3A5D3D", linetype = "dashed")+
labs(y = "Variables",#labels
x = "Importance (Scaled)")+
ggtitle("Random Forest Variable Importance")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 24, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 16),
axis.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "all_imp.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(aes(y = fit_all_frac$finalModel$y, x = fit_all_frac$finalModel$predicted), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", aes(y= fit_all_frac$finalModel$y, x = fit_all_frac$finalModel$predicted), color = "#3A5D3D")+#smooth line data
scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Elevation, Fractional Cover, and Indices")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 24, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 16),
axis.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "all_plot_caret.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
all_corr <- cor(y= fit_all_frac$finalModel$y, x = fit_all_frac$finalModel$predicted)
all_corr
## [1] 0.8453688
lm(fit_all_frac$finalModel$y ~ fit_all_frac$finalModel$predicted)
##
## Call:
## lm(formula = fit_all_frac$finalModel$y ~ fit_all_frac$finalModel$predicted)
##
## Coefficients:
## (Intercept) fit_all_frac$finalModel$predicted
## -0.03554 1.01109
max(fit_all_frac$finalModel$predicted)
## [1] 10.69074
min(fit_all_frac$finalModel$predicted)
## [1] 0.6010528
raster_all <- c(crsi_all, ndvi_all, vssi_all, mari_all, mllw_all, gv_frac_all, soil_frac_all, npv_frac_all)
all_pred_frac <- predict(raster_all, fit_all_frac, drop_dimensions = F, na.action = na.omit)
all_pred_crop <- all_pred_frac %>%
st_crop(y = boundary, crop = TRUE, as_points = FALSE) %>%
filter(date == "2022-02-24")
#write_stars(all_pred_crop, dsn = here("data", "predictions_nc_feb.tif"))
ggplot()+
geom_stars(data = all_pred_crop, aes(x = x, y = y, fill = prediction), na.action = na.omit) +
scale_fill_gradientn(colors = rev(cal_palette(name = "kelp1", n = 4, type = "continuous")))+
coord_sf(crs = 32116)+
labs(y = "Latitude",#labels
x = "Longitude",
fill = "Electrical
Conductivity
(mS/cm)")+
ggtitle("February 24, 2022")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"))
set.seed(1234)
fit_all_gv <- train(electro_cond_mS_per_cm ~ elevation + mari + vssi + GV_fraction,
data = soils_extract_2,
method = "rf",
trControl = trainControl(method = "cv", number = 10),
importance = TRUE)
# using caret
fit_all_gv_imp <- as.data.frame(fit_all_gv$finalModel$importance)
fit_all_gv_imp_scaled <- predict(preProcess(fit_all_gv_imp, method = c("range")), fit_all_gv_imp) %>%
rownames_to_column(var = "variable")
ggplot(fit_all_gv_imp_scaled) +
geom_point(aes(x = IncNodePurity, y = variable), color = "#E69512")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = IncNodePurity), color = "#E69512", size = 1.2)+
geom_point(aes(x = fit_all_gv_imp_scaled$`%IncMSE`, y = variable), color ="#3A5D3D")+
geom_segment(aes(y= variable, x = 0, yend = variable, xend = fit_all_gv_imp_scaled$`%IncMSE`), color = "#3A5D3D", linetype = "dashed")+
labs(y = "Variables",#labels
x = "Importance (Scaled)")+
ggtitle("Random Forest Variable Importance")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 24, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 16),
axis.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "all_imp.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
ggplot()+#ggplot object
geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
geom_point(aes(y = fit_all_gv$finalModel$y, x = fit_all_gv$finalModel$predicted), color = "#E69512")+#data to make points from
geom_smooth(method = "lm", aes(y= fit_all_gv$finalModel$y, x = fit_all_gv$finalModel$predicted), color = "#3A5D3D")+#smooth line data
scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
x = "Predicted Electrical conductivity (mS/cm)")+
ggtitle("Elevation, GV Fractional Cover, and Indices")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 24, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41", size = 16),
axis.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"))
#ggsave(filename = "all_plot_caret.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot
all_corr <- cor(y= fit_all_gv$finalModel$y, x = fit_all_gv$finalModel$predicted)
all_corr
## [1] 0.8405505
lm(fit_all_gv$finalModel$y ~ fit_all_gv$finalModel$predicted)
##
## Call:
## lm(formula = fit_all_gv$finalModel$y ~ fit_all_gv$finalModel$predicted)
##
## Coefficients:
## (Intercept) fit_all_gv$finalModel$predicted
## 0.0188 1.0022
max(fit_all_gv$finalModel$predicted)
## [1] 10.21696
min(fit_all_gv$finalModel$predicted)
## [1] 0.5007959
all_pred_gv <- predict(raster_all, fit_all_gv, drop_dimensions = F, na.action = na.omit)
plot(all_pred_gv)
all_pred_crop <- all_pred_frac %>%
st_crop(y = boundary, crop = TRUE, as_points = FALSE) %>%
filter(date == "2022-02-24")
#write_stars(all_pred_crop, dsn = here("data", "predictions_nc_feb.tif"))
ggplot()+
geom_stars(data = all_pred_crop, aes(x = x, y = y, fill = prediction), na.action = na.omit) +
scale_fill_gradientn(colors = rev(cal_palette(name = "kelp1", n = 4, type = "continuous")))+
coord_sf(crs = 32116)+
labs(y = "Latitude",#labels
x = "Longitude",
fill = "Electrical
Conductivity
(mS/cm)")+
ggtitle("February 24, 2022")+#plot title
theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
plot.background = element_rect("white"),
panel.background = element_rect("#faf7f2"),
panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
axis.text = element_text(color = "#5b4f41"),
axis.title = element_text(color = "#5b4f41"),
strip.background = element_rect("#f8f8f8"),
strip.text = element_text(color = "#5b4f41"),
axis.line = element_line(color = "#5b4f41"),
legend.title = element_text(color = "#5b4f41"),
legend.text = element_text(color = "#5b4f41"))